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I.  INTRODUCTION 

The  goals  of  this  project  are  to  insert  mature  diffractive  scattering  and  radiation  technology  into 
the  Xpatch®  computational  electromagnetics  (CEM)  prediction  codes  and  analysis  tools. 
Diffraction  plays  a  significant  role  in  the  radar  scattering  from  targets  consisting  of  smooth 
surfaces  or  joints,  and  numerous  applications  such  as  Reduced  Radar  Cross  Section  (RRCS) 
aircraft  design,  threat  signature  exploitation,  automatic  target  recognition  (ATR),  and  sensor 
fusion  (SF)  research  and  development  can  benefit  from  the  accuracy  improvements  associated 
with  diffractive  scattering  models.  The  diffraction  enhancements  developed  under  this  project 
shall  be  delivered  to  the  Air  Force  customer  as  well  as  offered  in  general  Xpatch®  releases  to  the 
United  States  Department  of  Defense  (DoD)  scattering  and  radiation  community  at  large  for 
multiple  wideband  applications. 

During  Year  1  (the  fourteen  months  from  July  2002  through  August  2003)  the  contractor 
integrated  a  robust  implementation  of  the  newly  developed  Air  Force  Office  of  Scientific 
Research  (AFOSR)  mathematical  formulations  for  higher-order  scattering  mechanisms  into 
Xpatch®  to  calculate  more  accurate  current  estimates  in  the  lit  as  well  as  the  shadow  region  of  the 
targets.  Specifically,  the  contractor  reformulated  truncated  wedge  incremental  length  diffraction 
coefficients  (ILDCs)  and  shadow  boundary  ILDCs  and  inserted  them  into  Xpatch®.  In  addition, 
the  supporting  Xpatch®  geometry  tools  were  modified  to  support  these  new  ILDC  capabilities. 
The  contractor  employed  appropriate  testing,  validation  criteria,  and  documentation  to  enable  the 
rapid  transition  of  this  software  technology  to  a  broad  range  of  United  States  DoD  applications. 

During  Year  2  (the  nine  months  from  September  2003  to  June  2004)  the  contractor  performed  a 
detailed  study  of  creeping  waves  to  supplement  the  shadow  boundary  ILDC  development  of 
Year  1.  This  study  was  motivated  by  the  role  of  creeping  waves  as  a  source  of  interaction 
between  installed  antenna  arrays  and  air  frames.  A  prototypical  creeping  wave  ray  tracer  was 
utilized  for  the  purposes  of  the  study.  Future  efforts  beyond  the  scope  of  this  project  will  focus 
on  development  of  a  robust  creeping  wave  ray  tracer. 
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II.  SHADOW  BOUNDARY  ILDC  DEVELOPMENT 


A.  Introduction 

The  contractor  has  implemented  the  shadow  boundary  ILDCs  in  both  their  simple  and  extended 
form  into  Xpatch®.  The  reformulation  of  the  detailed  expressions  of  the  new  shadow-boundary 
ILDCs  for  Xpatch®  has  been  the  primary  task  under  this  effort  with  close  collaboration  with  Dr. 
Arthur  Yaghjian.  The  underlying  implementation  of  the  ILDCs  in  Xpatch®  was  designed  such 
that  the  geometric  representation  of  the  underlying  curvature  and  shadow  boundary  of  the  target 
initial  graphics  exchange  specifications  (IGES)  CAD  file  are  used  to  calculate  all  the  necessary 
information  for  the  ILDCs,  without  requiring  intervention  by  the  Xpatch®  user. 

Without  the  added  capabilities  realized  with  shadow  boundary  ILDCs,  Xpatch®  forward 
scattering  degrades  for  smooth  structures.  Figure  1  compares  the  bistatic  scattering  from  a 
circular  cylinder  computed  by  Xpatch®  without  shadow  boundary  ILDCs  to  the  exact  solution. 
The  Xpatch®  solution  degrades  as  the  observation  angle  moves  from  the  lit  to  the  shadow  zone. 

To  illustrate  the  importance  of  the  shadow  boundary  ILDC  capability,  consider  the  imaging 
problem  illustrated  in  Fig.  2  of  a  tank  shadowed  by  a  hill’s  crest.  Currently  SAR  simulations  are 
line-of-sight  based.  Generic  scene  modeling  will  miss  the  tank  based  on  the  first-order 
assumption  that  no  incident  energy  reaches  the  tank.  However  the  true  radar  will  see  returns 
from  the  tank.  The  application  of  shadow  boundary  ILDCs  along  the  ridge  crest  permits 
Xpatch®  modeling  of  these  types  of  scenarios. 

B.  Technical  Overview 

Shadow  boundary  incremental  length  diffraction  coefficients  (SB-ILDCs)  are  not  a  complete 
solution  to  capturing  the  physics  embodied  in  the  diffraction  and  radiation  of  energy  in  the 
shallow  shadow  region.  However,  SB-ILDCs  do  represent  a  significant  improvement  over  the 
standard  physical  optics  (PO)  formulation  by  providing  a  far-field  correction  to  the  discontinuous 
PO  fields  along  the  shadow  boundary  (which  abruptly  step  from  finite  on  the  lit  side  to  zero  on 
the  shadowed  side).  In  effect  the  shadow  boundary  ILDCs  soften  the  shadow  and  improve  the 
bistatic  scattering  prediction  of  Xpatch®,  but  they  remain  inadequate  for  deep  shadow  region 
modeling  of  creeping  waves  traveling  long  distances  along  geodesic  paths. 
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Bistatic  Scattering  from  Circular  Cylinder 
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Fig.  1 .  A  6”  radius  cylinder  (top)  of  2”  height.  Bistatic  scattering  comparison  of  Xpatch  with  an  exact 
modal  solution  (bottom).  Elinc=Elobs=0°,  Azinc=0°,  18.8  GHz,  TMZ  polarization,  ka  =  60. 


Fig.  2.  Tank  shadowed  by  the  crest  of  a  hill.  A  geometrical  optics  or  physical  optics  approach  will  assume 
that  no  incident  energy  reaches  the  target  when  in  fact  radar  energy  is  diffracted  by  the  shadow  boundary 
formed  at  the  crest  of  the  hill  and  illuminates  the  target. 
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The  present  formulation  is  based  upon  the  ILDC  work.  [1].  The  approach  is  built  on  closed-form 
2-D  expressions  of  scattering  from  the  circular  cylinder  [2]  and  the  previous  ILDC  work  on  the 
prolate  spheroid  [3].  Non-uniform  (NU)  currents  in  the  vicinity  of  the  shadow  boundary  are  the 
difference  between  the  exact  currents  and  those  predicted  by  PO.  These  currents  are  significant 
in  the  immediate  vicinity  of  the  shadow  boundary.  The  shadow  boundary  ILDC  currents  are  the 
NU  currents  on  a  circular  cylinder  that  locally  conforms  to  the  surface  of  the  scatterer. 

The  geometric  parameters  needed  to  calculate  the  shadow  boundary  ILDCs  are  the  shadow 
boundary  curve,  the  vector  normal  to  the  shadow  boundary  surface  and  the  scatterer’ s  radius  of 
curvature  across  the  shadow  boundary.  This  short  list  of  parameters  can  be  efficiently  computed, 
providing  a  low  computational  cost  to  using  the  shadow  boundary  ELDCs. 

The  equations  for  the  shadow  boundary  ELDCs  are  given  below.  The  integral  of  dESB  along  the 
shadow  boundary  estimates  the  electric  field  due  to  the  NU  currents. 
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The  above  expressions  assume  a  harmonic  dependence  of  e  m  (trivial  transformations  need  to 


be  made  when  comparisons  to  Xpatch®  are  made  because  Xpatch®  assumes  a  e'c°‘  time- 
dependence).  The  wave  propagation  constant  in  free-space  is  given  by  k  =  a>/c ,  where  c  is  the 

speed  of  light  and  &>>  (K  The  subscript  ^  denotes  a  local  coordinate  system  (X/ ’^,z4  whose 
origin  is  a  point  SB  on  the  shadow  boundary.  The  positive  ■y‘  -axis  is  in  the  direction  of  the 
external  normal  UsB  to  the  surface  of  the  scatterer  at  Fsb  .  The  z,-axis  is  the  tangent  to  the 
shadow  boundary  curve  at  ^SB ,  with  a  positive  direction  assigned  to  the  Z/  -axis  by  requiring  that 
k  x,  k  (nSB  x  z, )  >  0  -phe  spherical  polar  angle  &01  is  the  angle  between  -  kmc  ancj  ^1 


The  angles  1  and  are  the  spherical  angles  of  the  direction  to  the  observation  point;  and 

are  the  corresponding  unit  spherical  angle  vectors;  r{  is  the  distance  from  PSB  to  the  observation 

point,  and  dz,  is  the  differential  length  along  the  shadow  boundary  at  PSB .  The  scalar  ac  is  the 
radius  of  the  perfectly  electric  conducting  (PEC)  circular  cylinder  whose  axis  is  parallel  to  the 
shadow  boundary.  The  functions  FzSB(kac,y/)  and  Gf  (kac , y/)  are  the  2-D  electric  and  magnetic 
far-fields  radiated  by  the  non-uniform  currents  at  the  top  shadow  boundary  of  the  cylinder 
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illuminated  by  a  TM  and  TE  plane  wave,  respectively,  of  unit  amplitude  normally  incident  in  the 
positive  x-direction.  The  coefficients  E™  and  Elf  are  the  complex  amplitudes  of  the  TM  and 
TE  components  of  the  incident  field  at  PSB ,  that  is 

E'’‘(PSB)  =  -E"0o,+E*y,  (13) 

where  E™  and  Eff  are  given  by 

U.  \  -  ,  \  (14) 

=  (k-‘xhst)-E-‘{PSB) 

and 

Eo,  =™sb'E  (Esb)  (15) 

and  0O!  =  hSB  x  kmc  and  y,  =  hSB .  The  summations  over  the  creeping  waves  indexed  by  s,  which 
continue  to  infinity,  have  been  truncated  to  S^.  This  implies  that  creeping  waves  are  not 
followed  or  modeled  deep  within  the  shadowed  region.  The  circular  cylinder  coordinates 
centered  on  this  cylinder  are  denoted  by  (p,y/)  (see  Fig.  3). 

Results  presented  in  the  next  section  of  the  report  demonstrate  the  effectiveness  of  this 
formulation.  In  summary,  the  main  assumptions  made  in  deriving  these  formulas  are  as  follows: 

1)  Plane- wave  incident  field 

2)  Constant  radius  of  curvature  in  the  longitudinal  direction  at  the  shadow  boundary 

3)  Infinite  radius  of  curvature  in  the  transverse  direction  at  the  shadow  boundary 

4)  Creeping  wave  truncation 

5)  Phase  correction  added  based  on  empirical  data 

C.  Shadow  Boundary  Search  Algorithm 

The  shadow  boundary  search  algorithm  calculation  effort  is  of  0(n),  which  is  similar  to  PO.  It 
requires  exhaustive  interrogation  of  the  surface  geometry  and  must  be  recalculated  for  every  new 
direction  of  incident  illumination.  The  calculation  is  performed  for  curved  surfaces  only.  The 
algorithm  first  breaks  the  surfaces  into  sufficiently  small  (typically  0.1  X)  patches.  The  shadow 
boundary  search  algorithm  then  computes  the  quantity  h  •  kinc  =  0  at  the  four  comers  of  the  patch. 

Figure  4  shows  a  surface  patch,  with  its  comers  labeled  with  a  blue  dot.  Finally,  the  algorithm 
interpolates  zero-crossings  to  estimate  the  shadow  contour.  This  procedure  locates  the  shadow 
boundary  contour,  where  h  ■  kjnc  =  0 .  A  sophisticated  root  finder  is  not  necessary  for  this 

calculation,  and  the  piecewise  linear  search  approach  converges  towards  the  true  shadow 
boundary  contour  as  the  discretization  size  is  made  finer. 
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Shadow  boundary  of 
a  circular  cylinder 


Fig.  3.  Shadow  boundary  quantities  relating  to  the  illumination  of  a  circular  cylinder.  The  local  coordinate 
system  is  defined  such  that  the  axis  of  the  cylinder  is  in  the  z-direction. 


Fig.  4.  Illustration  of  the  potential  illumination  scenarios  affecting  a  small  cell  of  a  curved  patch. 
The  shadow  boundary  is  approximated  by  a  line  segment  within  each  affected  cell. 


D.  Results  and  Discussion 

This  section  presents  computational  results  illustrating  the  efficacy  of  the  shadow  boundary 
ILDCs.  The  first  result  is  the  bistatic  (VV-pol)  scattering  from  a  circular  cylinder  illuminated 
perpendicular  to  its  axis.  The  boundary  search  algorithm  places  two  segments  on  the  cylinder’s 
shadow  boundary  as  shown  in  Fig.  5.  The  results  show  that  Xpatch®  with  the  shadow  boundary 
correction  is  a  significant  improvement  over  the  PO  solution  and  is  nearly  identical  to  the  exact 
solution. 

The  second  example  is  a  six-inch  radius  sphere.  Here,  the  shadow  boundary  is  modeled  by 
segments  located  along  the  equator  of  the  sphere.  Figure  6  identifies  the  segmented  shadow 
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boundary  on  the  sphere.  The  neighboring  graph  illustrates  the  accuracy  gained  by  including  the 
SB-ILDCs  in  the  computation.  The  Xpatch®  with  SB-ILDC  result  is  nearly  identical  to  the  exact 
solution  while  the  Xpatch®-PO  result  exhibits  significant  deviations. 

The  next  example  investigates  the  SB-ILDC  contribution  to  the  scattering  solution.  The  graphs 
shown  in  Fig.  7  compare  the  exact  minus  the  PO  solution  to  the  SB-ILDC  (only)  calculated 
scattering.  The  agreement  is  good  to  about  45  dB  down.  The  discrepancies  can  be  traced  to  the 
20  cells  per  wavelength  PO  approximation. 

The  SB-ILDC  plus  PO  solution  maintains  good  accuracy  compared  to  the  standard  Xpatch® 
solution  as  the  scatterer  size  decreases.  Figure  8  illustrates  the  scattering  from  a  sphere  as  its  size 
decreases  from  ka  =  20  to  ka  =  6.  Figure  9  shows  the  scattering  from  an  elliptical  cylinder.  This 
geometry  has  a  non-constant  radius  of  curvature,  suggesting  a  more  complex  field  variation  along 
the  shadow  boundary.  Figure  10  shows  the  scattering  from  a  spheroid  illuminated  by  a  plane 
wave  with  at  Az  -  -45°  incidence.  Internally,  Xpatch®  employs  Euler’s  formula  to  calculate  the 
local  radius  of  curvature  along  each  subdivision  of  the  shadow  boundary  contour. 
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Az  (deg) 

Fig.  5.  The  6”  radius  cylinder  (top)  of  Fig.  1  with  the  ILDC  segments  shown, 
comparison  of  Xpatch  with  shadow  boundary  correction  versus  an  exact  modal 
Elinc ~El0\yS — 0° ,  Aznc=0°,  18.8  GHz,  TMZ  polarization,  ka  =  60. 


Bistatic  scattering 
solution  (bottom). 


11 


Contract  No:  F49620-02-C-0038 


Horizontal  Polarization 


Fig.  6.  A  6”  radius  sphere  (top)  with  the  shadow  boundary  segments  shown.  Bistatic  scattering  comparison 
of  Xpatch®  with  shadow  boundary  correction  versus  an  exact  modal  solution  (bottom).  Elinc=Elobs=0°, 
Azmc=0°,  18.8  GHz,  TMZ polarization,  ka  ~  60. 
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Horizontal  Polarization 
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Fig.  7.  The  isolated  shadow  boundary  ILDC  contributions  (physical  optics  component  removed) 
for  the  sphere  of  Fig.  6. 
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Fig.  10.  Bistatic  scattering  from  the  prolate  spheroid  for  Az  =  -45°  incidence  with  ka  =  50  and  kb  =  25. 
The  shadow  boundary  determined  by  Xpatch  is  shown  in  light  green  and  agrees  with  the  closed  form 
solution  [3]  to  the  shadow  boundary  contour. 


E.  Summary 

The  implementation  of  shadow  boundary  ILDCs  in  Xpatch®  has  been  demonstrated  by  the 
improvement  of  scattering  predictions  for  large  cylinders  and  spheres  with  a  wide  range  of  sizes. 
Scatterers  with  changing  radii  of  curvature  were  successfully  modeled,  and  adequate  accuracy 
was  realized  even  for  rapidly  changing  radii  of  curvature.  The  implemented  shadow  boundary 
search  algorithm  was  verified  against  trivially  defined  shadow  boundaries  for  the  cylinder  and 
sphere,  as  well  as  for  more  complex  shadow  boundaries,  such  as  the  prolate  spheroid  at  off-axis 
incidence  angles. 
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III.  TRUNCATED  WEDGE  ILDC  DEVELOPMENT 
A.  Introduction 

The  contractor  has  reformulated  the  truncated  wedge  incremental  length  diffraction  coefficients 
(TW-ILDCs)  and  has  inserted  them  into  both  the  time-  and  frequency-domain  versions  of 
Xpatch®  in  a  manner  that  will  facilitate  their  use  for  complex  targets.  User-friendly  tools  for  the 
truncated  wedge  ILDCs,  (e.g.,  a  tool  to  extract  the  "other"  side  of  the  wedge)  as  well  as  a  robust 
and  accurate  electromagnetic  (EM)  implementation,  have  been  implemented  under  this  task. 
Validation  and  testing  of  the  truncated  wedge  ELDCs  have  been  accomplished  with  complex 
targets  represented  by  curved  surface  IGES  CAD  formats. 

The  sole  use  of  PO  to  calculate  the  scattering  from  various  PEC  bodies  results  in  significant 
inaccuracies  for  bodies  possessing  edges.  Inclusion  of  the  fringe  wave  (FW)  field,  a  field  that 
represents  diffraction  effects  due  the  edge,  improves  the  accuracy  of  such  calculations.  The  FW 
field  is  calculated  from  an  integral  along  incremental  strips  located  on  a  canonical  wedge  or  half¬ 
plane  [4-7],  The  standard  Xpatch®  code  calculates  these  wedge  ILDCs.  Symbolically  the  FW 
field  is  expressed  as  E^W  =  eGTD  _  gPO  for  the  semi-infinite  wedge.  The  fields  are  expressed 
as  a  semi-infinite  integral  along  the  projection  of  the  Keller  cone  on  the  wedge  face.  This 
integral  can  be  expressed  as  a  line  integral  along  a  differential  segment,  and  the  resulting 
contribution  is  referred  to  interchangeably  as  an  equivalent  edge  current  (EEC)  or  incremental 
length  diffraction  coefficient  (ILDC). 

The  use  of  truncated  wedge  ILDCs  is  motivated  by  the  following  problem  (see  Fig.  11).  The 
classical  wedge  ILDCs  formulated  in  [4-7]  possess  non-removable  “Ufimtsev”  singularities  when 
kinc  approaches  grazing  and  kobs  is  near  kinc-  The  classical  wedge  ILDCs  also  exhibit 
discontinuities  in  the  fringe  wave  field  due  to  the  use  of  non-physical  semi-infinite  current  strips. 


Fig.  1 1 .  Fringe  wave  currents  emanating  from  the  left  edge  of  a  four-sided  patch.  The  fringe  wave  currents 
are  strong  for  incident  fields  perpendicularly  polarized  relative  to  the  edge  while  fringe  wave  currents  due 
to  parallel  incidence  are  negligible. 

B.  Technical  Overview 

As  stated  above,  the  classical  wedge  ILDCs  exhibit  infinities.  Figure  12  shows  the  scattering 
from  an  equilateral  triangle  with  a  two-meter  side  length  and  a  one-meter  cylinder  height.  The 
scatterer  has  nine  diffracting  edges.  The  yellow  edge  is  nearly  at  grazing  incidence,  primarily 
along  the  skew  direction.  Surface  currents  in  the  direction  of  the  Keller  cone  are  non-decaying 
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[8].  A  singularity  exists  due  to  a  non-converging  integral.  ( kobs  ->  0  or  x  in  the  edge  local 
coordinate  system.) 


HH-pol  Bistatic  Scattering 


-40  -t . . — — . r ■ . . . . . — t . . J 
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Fig.  12.  Infinities  and  discontinuities  in  the  radar  cross-section  (RCS)  arising  from  the  single  orange  wedge 
ILDC  (top)  as  calculated  using  the  classical  ILDC  formulation.  The  physical  optics  component  of  the 
scattered  field  is  not  included  in  order  to  emphasize  the  behavior  due  to  wedge  diffraction. 


The  wedge  ILDC  formulation  implemented  in  Xpatch®  differs  from  the  rigorously  derived  wedge 
ILDCs  in  that  the  Keller-Ufimtsev  physical  theory  of  diffraction  (PTD)  coefficients  with  singular 
terms  are  omitted.  While  the  infinities  of  Fig.  12  are  avoided,  this  omission  of  terms  results  in 
numerous  discontinuities  in  the  computed  RCS  magnitude  and  phase.  Figure  13  shows  the 
discontinuities  in  the  predicted  RCS  of  a  triangular  cylinder. 
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Fig.  13.  Discontinuities  in  the  RCS  arising  from  the  three  orange  wedge  ILDCs  (top)  as  calculated  by 
Xpatch®  using  modified  form  of  classical  ILDC  formulation.  The  physical  optics  component  of  the 
scattered  field  is  not  included  in  order  to  emphasize  the  behavior  due  to  wedge  diffraction. 


C.  The  Truncated  Wedge  Algorithm 

In  order  to  improve  upon  the  wedge  ILDCs  used  in  the  standard  Xpatch®  code,  truncated  wedge 
ILDCs  are  employed.  Developed  by  Johansen  [8],  these  TW-ILDCs  are  well-behaved  for  all 
directions  of  incidence  and  observation  and  yield  a  finite  value  when  the  strip  length  goes  to 
zero. 

The  implementation  of  the  TW-ILDCs  is  a  follows.  First  each  edge  is  associated  with  its 
adjoining  faces.  As  shown  in  Fig.  14,  face  one  has  four  edges;  face  two  has  three  edges;  face 
three  has  four  edges,  and  so  forth.  Then,  differential  strips  are  traced  along  each  face,  and 
integrations  of  the  fringe  wave  contributions  along  the  finite  length  strips  are  calculated.  These 
finite  integrations  are  guaranteed  to  converge.  The  fringe  wave  current,  JFW  =  fTD  -  f°  is 
restricted  to  L  and  accounts  for  the  transition  region  (small  electrical  distances  from  the  edge) 
and  remains  valid  for  large  L  extending  well  beyond  the  transition  region. 

The  face  association  criteria  are:  shared  end  points,  at  least  one  shared  normal,  and  formation  of 
a  closed  polygon.  The  criteria  are  calculated  automatically  by  Xpatch®,  and  users  can  intervene 
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via  the  ModelMan®  CAD  tool  if  special  cases  such  as  interior  edges  (not  modeled  by  TW- 
ILDCs)  form  part  of  a  wedge  face. 

The  effort  required  to  calculate  the  TW-ILDCs  compared  to  classical  wedge  ILDCs  is  increased 
due  to  presence  of  Fresnel  integrals.  Also,  the  truncated  segments  vary  with  incidence  angle  and 
position,  requiring  new  calculations  for  each  angle  of  incidence  and  sample  point  along  the  edge, 
as  displayed  in  Fig.  15. 


Fig.  14.  Face  and  edge  associations  for  the  triangular  wedge  of  Fig.  12.  Edges  which  touch  (or  nearly  touch 
to  within  an  acceptable  tolerance)  and  possess  matching  normals  on  at  least  one  of  their  sides  are  considered 
to  be  associated.  Shown  on  the  right  in  orange  are  the  two  fringe  wave  path  segments  emanating  from  a  hit 
point  on  the  center  edge  and  traveling  until  they  reach  an  associated  edge. 


Fig.  15.  Illustration  of  how  varying  angles  of  incidence  result  in  different  current  strip  path  directions  and 
path  lengths.  Truncated  wedge  current  strips  are  a  function  of  hit  point,  incidence  angle,  arrangement  of 
neighboring  associated  edges,  as  well  as  frequency  and  wedge  angle. 


D.  Results  and  Discussion 

To  illustrate  the  efficacy  of  the  new  TW-ILDCs  the  scattering  from  a  triangular  cylinder  has  been 
calculated.  For  this  example  the  incident  angle  is  normal  to  the  edges  (see  Fig.  16),  with  three 
edges  illuminated.  A  comparison  is  made  between  the  old  ILDC  method,  the  truncated  wedge 
ILDC  method,  and  the  2-D  moment  method  solution  with  PO  subtracted  out.  Note  that  the 
truncated  wedge  solution  is  continuous  across  the  aspect  of  illumination.  The  non-physical 
infinities  associated  with  the  old  wedge  ILDCs  have  been  eliminated. 
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Fig.  16.  Comparison  of  the  new  truncated  wedge  ILDCs  for  the  wedge  configuration  of  Fig.  12.  The 
physical  optics  component  of  the  scattered  field  is  not  included  in  order  to  emphasize  the  behavior  due  to 
wedge  diffraction. 

Figure  17  illustrates  the  old  ILDC  and  new  TW-ILDC  scattering  from  a  plate  at  skew  incidence. 
This  example  is  a  square  PEC  plate  with  a  side  length  of  0.5  meters  (3  A.),  illuminated  at  a 
frequency  of  1 .799  GHz.  It  has  four  knife  edges  connected  to  two  faces.  The  incident  direction 
angles  are  elevation,  45  degrees  and  azimuth,  -45  degrees.  The  observation  azimuth  angle  is 
fixed  at  45  degrees,  with  the  elevation  ranging  from  0  to  180  degrees. 

A  final  example  is  the  “pencil”  geometry  of  Fig.  18.  This  case  was  originally  developed  at  the 
University  of  Illinois  during  studies  of  the  Fast  Illinois  Solver  Code  (FISC)  in  June  2003.  A 
circular  cylinder  3  m  in  length  and  0.1  m  radius  is  capped  at  one  end  with  a  circular  plate,  and  the 
other  end  terminates  in  a  circular  cone  of  length  0.167  m.  The  frequency  of  interest  is  3  GHz, 
and  the  incident  vector  is  directly  broadside  to  the  cylinder  (perpendicular  to  the  axis  of  the 
cylinder). 

The  pencil  of  Fig.  18  possesses  a  shadow  boundary  which  exists  along  the  cylinder  and  the  cone. 
The  principal  scattering  due  to  shadow  boundary  diffraction  is  expected  to  arise  from  the  cylinder 


21 


Contract  No:  F49620-02-C-0038 


because  the  shadow  boundary  along  the  cone  is  only  a  few  tenths  of  a  wavelength  in  length  and  is 
represented  by  segments  corresponding  to  diminishing  radii  of  curvature.  A  single  truncated 
wedge  face  is  placed  along  the  flat  circular  cap  at  the  other  end  of  the  pencil  geometry.  Shown  in 
Fig.  19  is  a  comparison  of  the  physical  optics  and  full  wave  results  for  the  bistatic  scattering  from 
the  pencil  geometry.  When  the  truncated  wedge  contributions  are  included,  the  agreement  with 
the  full  wave  solution  is  much  improved  (see  Fig.  20).  Finally,  when  the  shadow  boundary  ILDC 
contributions  are  included  in  addition  to  the  truncated  wedge  contributions,  further  agreement  is 
achieved  (see  Fig.  21).  The  improvements  due  to  the  shadow  boundary  ILDCs  will  be  strongest 
in  the  forward  scattering  direction.  Shown  in  Fig.  22  is  a  magnified  view  of  the  main  lobe  of  the 
pencil  scattering  in  which  the  shadow  boundary  contribution  is  necessary  to  correct  the 
magnitude  of  the  scattering  by  approximately  1  dB.  Neither  physical  optics  alone  nor  physical 
optics  with  the  truncated  wedge  contributions  are  capable  of  correctly  predicting  the  forward 
scattering  from  the  pencil  geometry. 


VV-pol 


Fig.  17.  A  square  plate  0.5  m  on  a  side.  Bistatic  scattering  comparison  of  Xpatch®  between  old  and  new 
ILDCs.  Elinc=Elobs=45°,Azi"c=-45°,  1.799  GHz. 


22 


RCS  (dBsm)  RCS  (dBsm) 


Contract  No:  F49620-02-C-0038 


Fig.  18.  The  pencil  geometry  composed  of  a  cone  and  cylinder  and  flat  circular  end  cap  (hidden). 
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Fig.  19.  Bistatic  scattering  from  the  pencil:  Full  wave  vs.  physical  optics. 


Fig.  20.  Bistatic  scattering  from  the  pencil  with  truncated  wedge  contributions  included. 
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Fig.  2 1 .  Bistatic  scattering  from  the  pencil  with  shadow  boundary  and  truncated  wedge  contributions  included. 
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Fig.  22.  Shadow  boundary  influence  in  forward  scattering  region  for  the  pencil  geometry.  The  shadow 
boundary  contribution  provides  a  significant  correction  at  the  purely  forward  scattering  direction 
(Az  =2  7(f) 

E.  Summary 

The  TW-ILDCs  significantly  improve  the  calculation  of  the  scattering  from  bodies  with  edges. 
Singularities  in  the  current  integration  associated  with  the  previous  wedge  ILDC  formulation 
have  been  removed,  and  discontinuities  that  were  present  in  the  calculated  scattered  field  are  no 
longer  present.  The  method  partially  takes  into  account  second  order  diffraction  effects.  The 
TW-ILDCs  derivation  includes  cross-polarization  terms  and  improves  reciprocity. 
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IV.  GEOMETRY  TOOL  DEVELOPMENT 

A.  Introduction 

The  existing  edge  and  ILDC  capabilities  in  ModelMan®  have  been  extended  to  handle  wedge 
face  definitions  appropriate  for  truncated  wedge  ILDC  analyses  in  Xpatch®.  A  simple,  visual 
interface  has  been  developed  for  automatically  creating,  modifying,  and  deleting  a  variety  of 
indexed  wedge  face  definitions.  Wedge  face  definitions  can  be  defined  via  three  methods  (i) 
direct  calculations  by  built-in  routines,  (ii)  import  from  existing  edge  files,  or  (iii)  manual  entry 
by  the  user.  Wedge  face  definitions  are  included  by  default  when  ModelMan®  saves  an  edge 
file. 

B.  Graphical  User  Interface  Features 

Shown  in  Fig.  23  is  the  ModelMan®  Wedge  Faces  tab.  This  basic  operations  available  through 
this  tab  are  similar  to  those  of  the  classical  edges  but  are  extended  to  include  settings  unique  to 
truncated  wedge  faces.  The  functions  of  this  tab  are  tabulated  below: 

Visibility  -  view  individual  wedge  faces  by  index,  usually  for  manual  customization  or 
deletion  (see  below),  or  all  wedge  faces  that  would  be  exported. 

Show  -  in  the  viewer,  show  representations  of  any  combination  of  the  wedge  faces,  edges  not 
associated  with  any  wedge  face,  and  faceted  geometry. 

Settings  -  sliders  for  automatic  generation  of  wedge  faces.  Wedge  angle  defines  minimum 
angle  between  adjacent  faces,  and  planar  tolerance  defines  the  maximum  angle  between  adjacent 
facets  on  the  same  face.  Generation  of  the  wedge  faces  is  dynamically  shown  in  the  viewer  as 
the  sliders  are  moved.  The  sum  of  the  two  angles  is  not  allowed  to  exceed  180. 

Selection  Options  -  allows  manual  customization  of  wedge  faces.  Once  automatic 
generation  of  wedge  faces  is  done  per  above,  individual  edges  can  be  removed  from  or  added  to 
the  current  face  by  picking  on  them  in  the  viewer.  Edges  are  not  allowed  to  belong  to  more  than 
two  faces. 

Buttons  -  Add  New  Face  -  adds  a  new  indexed  face  that  the  user  can  manually  populate  with 
edges.  Delete  Current  Face  -  deletes  the  current  wedge  face  definition,  leaving  the  edges  to  be 
associated  with  other  wedge  faces  if  desired. 
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Fig.  23.  The  ModelMan®  Wedge  Faces  tab.  Controls  for  visibility,  filtering  tolerances,  and  selection  of 
wedges  are  provided  for  through  this  interface. 


A  sample  geometry  containing  multiple  edges  is  shown  in  Fig.  24.  The  wedge  face  tab  is 
employed  to  identify  candidate  wedge  faces.  In  the  case  of  automatic  wedge  face  identification, 
the  user  has  the  option  to  disallow  any  wedge  faces  which  are  deemed  unacceptable  (such  as 
cases  arising  from  gaps  and  joints  which  give  rise  to  fictitious  edges).  Shown  in  Fig.  25  are  the 
outlines  of  the  wedge  faces  corresponding  to  the  original  geometry  of  Fig.  24. 
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V.  SHADOW  BOUNDARY  AND  CREEPING  WAVE  STUDY 
A.  Introduction 

One  of  the  key  assumptions  in  the  development  of  the  shadow  boundary  incremental  length 
diffraction  coefficients  (SB-ILDCs)  is  the  truncation  of  the  creeping  wave  in  the  deep  shadow 
region.  As  identified  previously,  the  SB-ILDCs  are  not  a  complete  solution  to  capturing  the 
physics  associated  with  the  shadow  region,  but  they  nonetheless  represent  a  significant 
improvement  over  the  physical  optics  formulation. 

In  this  component  of  the  project,  we  studied  practical  techniques  for  overcoming  this  deficiency 
of  SB-ILDCs.  In  particular,  techniques  were  investigated  for  extending  creeping  wave  rays  from 
the  shadow  boundary  (SB)  through  surface  ray  tracing  on  arbitrary  shapes  and  computing 
associated  surface  currents  to  augment  the  first-bounce  PO  solution  in  Xpatch®.  This  correction 
is  particularly  important  for  the  prediction  of  time-domain  radar  signatures  and  SAR  images  on 
smooth  targets,  since  energy  can  propagate  around  the  target  into  the  deep  shadow  region  and 
back  into  the  lit  region,  generating  delayed  response  features  that  are  missing  in  the  SB-TT.DC 
formulation. 

In  executing  this  work,  a  research  code  developed  at  the  University  of  Kentucky  was  employed  to 
accelerate  progress  on  the  studies.  This  code  also  implements  a  novel  solution  to  the 
identification  of  SBs  on  smooth  targets.  As  part  of  this  work,  we  studied  the  specifics  and 
performance  of  this  SB-finding  algorithm  and  have  compared  it  to  the  practical  algorithm  based 
on  ray  tracing  developed  for  the  incorporation  of  SB-ILDCs  into  Xpatch®. 

In  electromagnetic  scattering  problems,  the  PO  method  is  a  quick  way  of  determining  the  surface 
currents  on  the  scattering  body  without  the  need  to  rigorously  solve  Maxwell’s  equations.  This 
method  is  generally  suitable  when  the  target  and  key  target  features  are  electrically  large.  In 
practice,  target  features  on  the  order  of  a  wavelength  still  do  very  well  with  the  PO 
approximation,  and  larger  features  perform  even  better. 

One  of  the  well-known  anomalies  of  PO  is  that  it  only  predicts  reasonable  surface  currents  in  the 
lit  region  of  the  target,  while  non-physically  yielding  zero  surface  currents  on  surfaces  blocked 
from  direct  illumination.  In  other  words,  the  PO  surface  currents  abruptly  end  at  the  shadow 
boundary  (SB)  of  the  target  (with  respect  to  the  source),  even  if  the  target  surface  is  smooth  in  the 
vicinity  of  the  SB.  Among  other  problems,  this  creates  a  non-physical  scattering  feature,  which 
can  be  observed  and  localized  to  the  SB  in  the  time-domain  scattering  response.  In  reality,  the 
surface  currents  should  smoothly  transition  in  crossing  the  SB  from  the  lit  side  to  the  dark  side. 
Once  in  the  dark  side,  the  surface  currents  should  rapidly  decrease  as  they  shed  energy  to  free 
space  since  they  are  not  reinforced  through  direct  illumination. 

The  effect  is  illustrated  in  Fig.  26,  which  gives  the  time-domain  monostatic  backscattering  of  a 
PEC  sphere  (radius  a  =  0.0477  m)  based  on  the  exact  solution  to  Maxwell’s  equations  (Mie 
Series  [2])  and  the  PO  approximation.  In  the  figure,  the  time-axis  has  been  replaced  with  down 
range  distance  using  the  usual  technique  of  mapping  the  monostatic  response  time  delays  to  one¬ 
dimensional  positions  on  the  target  along  the  incident  /  back-scatter  axis.  The  initial  response  at 


28 


Contract  No:  F49620-02-C-0038 


z  =  -  a  =-0.048  m  corresponds  to  the  leading  edge  of  the  sphere,  and  here  PO  and  the  exact 
solution  are  in  excellent  agreement.  Agreement  on  this  critical  feature  in  the  time  response 
summarizes  in  a  special  case  the  popularity  of  the  PO  approximation  for  general  scattering 
problems  where  exact  solutions  are  not  available.  However,  beyond  this  leading  edge  response, 
the  exact  and  PO  solutions  markedly  diverge,  albeit  at  somewhat  lower  levels  than  the  initial 
response.  First,  the  PO  response  at  down  range  position  z  =  0  is  the  spurious  scattering  feature 
introduced  by  the  abrupt  and  non-physical  termination  of  surface  currents  at  the  sphere  SB. 
Second,  the  exact  solution  contains  a  delayed  response  at  z  =  m/2  =  0.075  m  that  is  completely 
absent  from  the  PO  solution.  This  is  caused  by  surface  tangential  creeping  waves  that  originate 
from  the  shadow  boundary  and  travel  around  the  sphere  along  its  meridians. 

The  creeping  waves  constantly  shed  energy  to  space  in  their  direction  of  travel,  and  when  they 
again  reach  the  SB  on  the  opposite  side  of  the  sphere  from  where  they  started,  their  radiation  is 
directly  in  the  back-scattering  (monostatic)  direction.  This  also  means  that  creeping  waves  can 
play  a  significant  role  in  bistatic  scattering  at  observation  angles  where  the  specular  (PO) 
scattering  mechanisms  are  weak. 


Range  Profile  of  PEC  Sphere  (a  -  0.0477  m) 


Fig.  26.  Monostatic  time-domain  response  (range  profile)  of  PEC  sphere  (a  =  0.0477  m):  physical  optics  and  Mie 
series  solutions.  Note  the  erroneous  late  time  responses  in  the  red  PO  curve  due  to  (a)  discontinuous  PO  fields  at  the 
shadow  boundary  and  (b)  omission  of  creeping  waves  traveling  in  the  deep  shadow. 


The  effect  of  the  discrepancy  between  PO  and  the  exact  solution  for  the  PEC  sphere  can  also  be 
observed  in  the  frequency  response,  as  shown  in  Fig.  27.  At  higher  frequencies  (koa  »  A),  PO 
and  the  exact  solution  converge  to  the  same  RCS  level,  and  this  is  to  be  expected  from  an 
asymptotic  solution  like  PO.  However,  the  variation  of  the  response  with  frequency  as  modeled 
by  PO  is  incorrect,  and  this  can  be  attributed  to  both  poor  approximation  of  surface  currents  in 
the  vicinity  of  the  SB  and  the  complete  absence  of  deep  shadow  surface  currents  associated  with 
the  creeping  wave  mechanism 
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B.  Technical  Overview 

The  core  methodology  of  Xpatch®  is  the  shooting-and-bouncing  ray  (SBR)  method,  which 
extends  the  lst-bounce  PO  solution  to  multi -bounce  through  geometrical  optics  (GO)  ray  tracing 
and  successive  application  of  PO  to  approximate  surface  currents  beyond  the  lst-bounce 
according  to  the  incident  and  reflected  GO  fields. 


Frequency  Response  of  PEC  Sphere 


Fig.  27.  Monostatic  frequency  response  of  PEC  sphere:  physical  optics  and  Mie  series  solutions.  Note  that  shadow 
boundary  ILDCs  will  provide  improvements  to  the  PO  for  ka>  6  (see  Fig.  8). 

More  abstractly  speaking,  Xpatch®  seeks  to  predict  electromagnetic  scattering,  and  associated 
effects  of  specular  reflection  and  diffraction,  through  an  asymptotic  approximation  of  surface 
currents  on  complex  targets.  While  Xpatch®  employs  geometric  ray  tracing  to  approximate  the 
currents,  the  terminal  point  (or  escape  direction)  of  tracked  rays  is  not  considered  in  establishing 
the  dominant  directions  of  scattered  energy.  Rather,  it  relies  on  the  radiated  wave  interference  of 
the  approximated  surface  currents  to  reconstruct  the  spatial  and  angular  regions  of  strong 
scattering.  In  this  respect,  SBR  corresponds  to  the  actual  (physical)  process  of  scattering.  One  of 
the  benefits  of  this  approach  is  that  complex  and  arbitrary  scattering  geometries  can  be  handled 
readily  without  the  need  to  search  for  dominant,  geometrically  minimized  ray  paths  terminating 
at  or  near  scattering  angles  of  interest.  In  developing  a  technique  for  handling  the  deep  shadow 
region  in  Xpatch®,  we  have  sought  to  evaluate  an  approach  that  is  consistent  with  the  underlying 
surface  current  approximation  framework  of  Xpatch®.  To  that  end,  the  technical  approach  can  be 
summarized  as  follows: 

1)  Identify  the  SB  on  the  target. 

2)  Sample  the  SB  uniformly  at  sub-wavelength  intervals. 

3)  Launch  tangential  creeping  waves  from  the  SB  into  the  shadow  region. 

4)  Trace  the  creeping  waves  over  smooth  target  surfaces  to  form  geodesics  (i.e., 
paths  that  locally  minimize  path  length  as  they  extend). 

5)  Generate  surface  currents  associated  with  the  found  creeping  wave  paths. 
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6)  Radiate  these  surface  currents  to  space  and  combine  the  result  with  the  base 
Xpatch®  PO-SBR  solution. 


1.  Finding  the  Shadow  Boundary  (Alternate  Method) 

By  definition,  the  shadow  boundary  is  the  locus  of  illuminated  points  on  the  target  where 
n  •  kinc  =  0,  which  we  refer  to  has  the  incident-normal  dot  product.  The  research  code  employed 

in  these  studies  implements  an  algorithmic  Method  B  for  finding  the  shadow  boundary  that  is 
somewhat  different  from  the  algorithmic  Method  A  presented  earlier  in  the  SB-ELDC  component 
of  this  project. 

Both  Methods  A  and  B  are  intended  for  targets  described  by  curved  surfaces.  Without  going  into 
great  detail  since  it  has  already  been  presented,  Method  A  employs  an  exhaustive  search  by 
looping  through  the  curved  surfaces  comprising  the  target  and  subdividing  them  into  small 
patches,  nominally  0.1  X  on  a  side.  Each  small  patch  is  nearly  planar  and  has  four  sides  and  four 
comers.  The  incident-normal  dot  product  is  computed  for  each  of  the  four  comers.  Based  on  the 
signs  of  these  dot  products,  it  can  be  established  whether  the  SB  intersects  the  patch.  If  it  does, 
the  approximate  path  of  the  SB  across  the  patch  can  be  established  through  linear  interpolation. 

Method  B  is  based  on  the  observation  that  the  SB  forms  a  closed  curve  along  the  geometry  of  the 
surface.  Method  B  consists  of  two-steps: 

1)  Find  one  point  on  the  target  along  the  SB.  This  point  designates  the  SB  start 
point. 

2)  Extend  the  SB  curve  along  the  target  surface  in  small  incremental  steps.  Repeat 
step  2  and  stop  when  the  SB  closes  back  on  itself  at  the  start  point. 

To  find  an  SB  start  point,  all  the  curved  surfaces  comprising  the  target  are  subdivided,  with  the 
limited  objective  of  yielding  mildly  curved,  quad  patches  that  do  not  fold  back  on  themselves.  A 
4x4  subdivision  is  usually  adequate.  The  algorithm  then  loops  through  these  patches  until  one 
is  found  with  a  side  that  intersects  the  SB.  This  determination  is  made  using  the  same  technique 
as  in  Method  A  of  testing  the  incident-normal  dot  product  at  each  comer  of  the  patch.  A  sign 
change  in  the  dot  product  indicates  an  SB  intersection  with  the  side  bounded  by  the  two  comers. 
The  precise  location  of  the  SB  along  the  side  is  then  determined  by  uniformly  sampling  the  dot 
product  along  the  side  and  further  localizing  the  dot  product  zero-crossing  through  linear 
interpolation.  This  yields  an  SB  start  point  on  the  “start  patch,”  referred  to  as  SB -PI  (Shadow 
Boundary  Point  1).  The  sampling  and  interpolation  process  is  more  conveniently  implemented  in 
the  2-D  parameterized  space  of  the  patch  domain  (i.e.,  u,v  coordinates)  rather  than  directly  in  the 
3-D  target  CAD  space  of  patch  range  (i.e.,  x,y,z  Cartesian  coordinates). 

If  the  SB  intersects  one  side  of  a  patch,  it  must  intersect  another  side  as  well  at  SB-P2,  as  shown 
in  Fig.  28.  This  second  side  is  identified  through  the  dot  product  sign  testing.  The  SB 
intersection  with  the  second  side  is  then  localized  to  SB-P2  as  before.  Having  established  SB -PI 
and  SB-P2  on  the  patch  sides,  the  SB  can  now  be  extend  through  the  patch  from  SB -PI  on  the 
first  side  and  toward  SB-P2  on  the  second  side.  In  the  implementation  utilized  during  this  study, 
the  SB  contour  was  approximated  as  cutting  a  linear  path  in  the  2-D  parameterized  domain  space 
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of  the  patch,  and  this  was  found  to  be  adequate  for  all  test  cases.  While  higher-order  (i.e.,  more 
accurate)  means  of  extending  the  SB  through  the  patch  can  be  developed,  the  same  result  can  be 
achieved  by  employing  finer  subdivisions  in  the  piecewise  linear  approach  presented  above. 
Regardless,  the  present  approach  yields  many  SB  points  over  the  patch,  forming  a  curved  1-D 
path  in  the  3-D  target  geometry  space.  These  patch  interior  SB  points  are  not  enumerated  in  the 
notation  here,  but  they  are  collected  consecutively  in  the  algorithm. 


Fig.  28.  Illustration  of  Method  B  for  finding  the  shadow  boundary.  The  result  of  Method  B  is  a  contiguous  shadow 
boundary  spanning  entire  patches  and  encircling  the  entire  geometry. 


Having  established  the  SB  contour  along  a  single  patch  in  the  target,  the  next  step  is  to  extend  the 
SB  to  the  neighboring  patch  that  shares  the  same  side  containing  SB-P2  on  the  first  patch.  This 
can  be  reliably  achieved  with  a  topologic  map  of  the  target  geometry  patches.  The  mapping 
establishes  connectivity  along  the  fours  sides  of  each  patch  with  its  neighboring  patches.  Details 
of  the  construction  of  the  topologic  map  will  be  presented  later.  Once  the  first  side  of  the 
neighbor  patch  is  identified,  SB-P2  is  carried  over  by  appropriate  conversion  of  the  parameter 
coordinates  (u,v)  from  the  first  patch  to  the  second  patch. 

At  this  point,  the  Method  B  algorithm  is  now  in  exactly  the  same  state  for  SB-P2  on  the  second 
patch  as  it  found  itself  after  having  identified  SB-P1  on  the  first  patch.  Hence  the  process  can  be 
repeated,  and  the  SB  extended  across  the  target  geometry  from  patch  to  patch.  The  process 
completes  when  the  extended  SB  closes  back  on  SB-P1. 

An  example  of  the  results  of  this  algorithm  are  shown  in  Fig.  29.  After  the  SB  has  been 
constructed  as  a  list  of  points,  it  can  be  uniformly  re-sampled  through  a  spline  fit  to  achieve  any 
desired  density  of  points.  This  is  useful  in  establishing  launch  points  for  creeping  wave  rays  of 
uniform  width,  as  discussed  later.  The  whole  process  proves  to  be  computationally  efficient.  For 
the  example  given  in  Fig.  29,  on  the  order  of  a  CPU  second  was  required  to  establish  the  SB. 
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Fig.  29.  Shadow  boundary  (colored  points)  found  by  the  alternate  Method  B  on  a  cylinder  with  hemisphere  end 
caps.  The  incident  angles  is  60°  off  the  cylinder  broadside.  The  incident  far  field  ray  trace  is  also  shown  to  illustrate 
the  consistency  of  the  extracted  SB  with  the  boundary  of  the  incident  ray  field  on  the  target. 


As  mentioned,  continuing  the  shadow  boundary  from  one  curved  surface  patch  to  the  next  is 
effected  through  a  topological  map.  This  map  can  be  built  once,  after  the  original  subdivision  of 
the  target  curved  surfaces  into  smaller  patches.  In  the  code  used  for  this  study,  the  topologic  map 
was  constructed  by  evaluating  the  3-D  position  of  the  comers  on  all  the  patches  and  finding 
matches  among  these  comers.  By  allowing  a  small  tolerance  in  matching  the  comers,  it  is  not 
necessary  for  the  sides  of  neighboring  patches  to  perfectly  align.  This  is  an  important 
consideration  since  target  models  based  on  curved  surfaces  often  suffer  from  this  problem  to 
some  degree. 

In  the  examples  we  considered,  all  of  the  surface  patches  aligned  with  neighbors  at  their  quad 
comers.  In  practical  target  CAD  models,  it  is  common  to  encounter  patches  that  share  a  portion 
of  their  side  without  sharing  a  comer  (i.e.,  the  comer  of  one  patch  resides  in  the  middle  of  a  side 
of  its  neighbor  patch).  In  this  case,  comer  matching  will  not  work  in  constmcting  the  topology 
map,  and  more  rigorous  search  techniques  are  required.  One  practical  possibility  is  to  use  a 
geometric  ray  tracer  to  systematically  interrogate  the  space  in  the  immediate  vicinity  around  the 
patch  sides.  This  process  will  construct  a  complete  list  of  all  neighboring  patch  relations  and  can 
then  used  to  expedite  the  construction  of  the  topology  map.  As  with  the  comer-matching 
approach,  the  more  rigorous  ray  tracing  approach  to  identifying  neighbor  patch  connectivity  can 
be  implemented  with  a  small  spatial  tolerance  to  allow  for  inadvertent  gaps  between  patches. 

2.  Qualitative  Comparison  of  Shadow  Boundary  Search  Method 

Having  described  two  methods  for  finding  the  SB,  we  now  compare  their  relative  merits. 
Method  A,  based  on  exhaustive  search,  has  the  benefit  of  inherent  robustness.  On  a  complex 
target,  it  will  find  all  the  shadow  boundaries.  In  contrast,  Method  B,  as  described,  will  only  find 
one  of  the  contiguous  SBs.  This  is  because  it  starts  by  finding  one  SB  point  and  then  extending 
the  SB  curve  until  it  closes  on  itself.  On  all  but  simple  convex  targets,  there  are  multiple  SB 
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contours.  This  does  not  eliminate  Method  B  and  instead  suggests  some  mechanism  must  be 
added  to  seed  the  multiple  SBs. 

Method  B,  based  on  seeding  an  SB  and  extending  it  to  closure,  has  the  advantage  of 
automatically  yielding  a  contiguous  shadow  boundary.  In  contrast,  Method  A  produces  disjoint 
samples  of  the  SB.  For  implementation  of  the  SB-ILDCs,  this  is  not  an  issue,  but  it  would  be  a 
problem  in  launching  creeping  wave  rays  in  a  spatially  controlled  fashion.  This  does  not 
eliminate  Method  A  from  the  creeping  wave  application,  but  it  means  some  mechanism  must  be 
added  by  which  the  disjoint  points  can  be  grouped  in  to  one  or  more  contiguous  SBs  to  facilitate 
uniform  re-sampling. 

Finally,  with  Method  A,  the  computational  effort  is  proportional  to  the  surface  area  of  the  target 
measured  in  wavelengths  squared.  In  contrast,  with  Method  B,  the  computational  effort  is 
proportional  to  a  combination  of  the  number  of  target  model  patches  and  the  length  of  the 
shadow  boundary  contours.  In  terms  of  computational  scaling,  Method  B  would  appear  to  be 
better.  However,  this  does  not  mean  that  Method  B  will  always  be  faster  than  Method  A.  It  just 
means  that  for  a  large  enough  problem  size,  Method  B  should  be  faster  than  Method  A.  The 
actual  performance  of  these  methods  has  not  been  tested  at  this  reporting,  so  no  conclusions 
beyond  scaling  properties  are  made  at  this  time. 

In  closing,  we  mentioned  earlier  that  Method  B  would  require  some  means  of  seeding  multiple 
SBs  on  a  practical  target.  One  possibility  is  to  construct  the  first  SB  and  then  begin  the  process 
all  over,  searching  through  the  patches  to  find  another  incident-normal  dot  product  zero-crossing. 
Instead  of  stopping  automatically  on  first  zero-crossing,  the  search  continues  until  a  patch  side 
zero-crossing  is  found  that  is  not  already  part  of  an  existing  SB  contour.  On  the  final  iteration, 
the  search  will  yield  no  zero-crossings  that  are  not  previously  covered,  indicating  that  all  SBs 
have  been  found.  An  implementation  of  this  approach  appears  practical,  but  whether  or  not  it 
would  be  reliable  has  not  been  tested. 

C.  Method  for  Tracing  Creeping  Waves 

As  described  previously,  the  motivation  for  this  component  of  the  project  is  to  evaluate  practical 
means  of  obtaining  surface  currents  in  the  deep  shadow  region  of  targets  by  including  a  creeping 
wave  mechanism  to  predict  the  deep  shadow  surface  currents.  This  requires  launching  creeping 
wave  rays  from  the  shadow  boundary  and  tracing  them  on  the  target  surface  into  the  shadow 
region,  possibly  traveling  far  enough  to  emerge  back  into  the  lit  region  on  some  other  part  of  the 
target.  The  intent  is  to  use  the  creeping  rays  paths  as  domains  for  creeping  wave  surface  currents, 
which  are  themselves  propagated  according  to  the  geometrical  properties  of  the  creeping  rays  and 
the  underlying  surface. 

The  launch  points  are  established  by  uniformly  sampling  the  target  SBs,  as  discussed  earlier.  The 
launch  direction  is  the  direction  of  the  incident  wave  at  the  SBs,  which  by  definition  is  tangent  to 
the  target  surface.  As  required  by  the  asymptotic  creeping  wave  physics,  to  continue  the  rays 
along  the  target  surface,  they  must  maintain  a  geodesic  path.  Any  1-D  arc  conformal  to  a  smooth 
2-D  surface  and  occupying  a  3-D  space  has  the  property  of  arc  length,  which  is  simply  the 
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distance  traversed  along  the  arc.  For  a  geodesic  path  between  two  end  points  with  arc  length  so, 
any  small,  surface-conformal  perturbation  of  the  arc  between  the  end  points  would  yield  a  new 
arc  length  so  +  &,  where  fis  >  0.  This  is  the  definition  of  a  geodesic  curve  [10,1 1]. 

The  concept  of  a  geodesic  arc  is  simple.  Differential  equations  for  tracing  one  along  an  arbitrary 
curved  surface  are  somewhat  more  complicated  but,  in  the  end,  manageable  in  a  numerical 
implementation.  For  a  curved  surface  patch,  the  3-D  displacement  vector,  r,  to  a  point  on  the 
surface  is  parameterized  in  2-D  u,v  coordinates:  r  -  r(u,v).  As  shown  in  Fig.  30,  the  numerical 
geodesic  tracing  problem  is  as  follows:  given  point  Po(uo,vo)  on  the  surface  and  initial  direction 
t0  tangent  to  the  surface,  advance  the  curve  to  unknown  point  Pn(un,vn)  by  taking  small  steps  Au, 

Av  (also  unknown)  while  maintaining  a  geodesic  path  at  all  points  along  the  way.  To  be  clear, 
although  the  geodesic  path  resides  in  3-D  target  geometry  space,  it  will  be  advanced  through 
operations  in  parameterized  2-D  u,v  space. 


Fig.  30.  A  geodesic  path  started  at  known  point  P0  in  known  direction  f0  and  numerically  advanced  through  small 
incremental  steps. 


Since  the  1-D  geodesic  arc  lies  on  the  2-D  surface,  any  point  along  the  arc  can  be  expressed  in 
u,v  coordinates,  but  u  and  v  are  not  independent  variables  along  the  1-D  arc.  The  approach  for 
forming  the  geodesic  arc  is  to  numerically  solve  the  differential  equations  relating  the  u  and  v. 
These  equations  are  given  by  [10,1 1]  as 


(17) 
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The  coefficients  J4 n  ,  I2  n,  etc.  (called  Christoffel  symbols)  relate  to  differential  properties  of  the 
2-D  surface  and  are  given  later.  These  values  can  be  determined  directly  from  the  algebraic 
expressions  for  r(u,  v)  for  the  surface  patch.  At  any  given  point,  only  one  of  these  two  equations 
needs  to  be  solved,  the  choice  depending  on  the  current  state  du/dv  and  stability  considerations. 


From  t0 ,  one  can  determine  duo/dvo,  as  will  be  explained  shortly.  Let  us  assume  that  duo/dvo  <  1 

for  now.  Also,  we  can  scale  Auq  and  Avo  to  correspond  to  the  desired  step  size  in  3-D  space,  Aso- 
This  yields  four  parameters  for  P0:  u0,  v0,  Au0,  Av0. 

2  2 

From  (16),  d  u/dv  at  Po:uo,vo,Auo,Avo  is  determined  through  direct  substitution.  Implementing 
the  2nd-order  Euler  method,  we  obtain 

M,  =  u0  +  Am0, 

Vi  =  v0  +  Av0 , 

Amj  =  Au0  +~y  Lo,v0  (Av0)2, 
dv 

Avj  =  Av0 . 

The  general  iteration  becomes 


uM=u,+Ault  (19a) 

v/+i=v,+Av,,  (19b) 

A um  =  Aui  +  TT  Uw  (Av«  )2 »  (19c) 

dv 

Av;+1  =  Av . .  (19d) 

Proceeding  in  this  way,  it  may  happen  that  Aut  >  Av{  (i.e.,  Aui/Avl  <  1)  for  some  i.  When  this 
occurs,  numeric  stability  is  enhanced  by  switching  to  (17)  and  iterating  as 

uM  =  ui  +  Am,.  ,  (20a) 

vl+i  =v,+Av,.,  (20b) 

A  ui+l  =  Am,.  ,  (20c) 

Avi+i  =  Av,-  +  -y-j  |UI>;  (Am,)2  .  (20d) 


(18a) 

(18b) 

(18c) 

(18d) 


Higher-order  numeric  iteration  solutions  to  the  differential  equations  of  (16)  and  (17)  may  also  be 
implemented  to  improve  accuracy  or  allow  a  larger  step  size  with  each  iteration,  but  the  2nd-order 
Euler  is  the  simplest. 
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The  original  duo/dvo  at  Po  can  be  obtained  from  t0  in  the  following  manner.  For  a  1-D  arc 

conformal  to  the  surface  and  parameterized  as  r  =  r(u,v)  =  r(u(t),v(t)),  the  tangent  velocity  vector 
is  given  by 

-  dr  _  du  _  dv 

t~~dt~Vu~dt+rv~di'  (21) 

Here,  we  use  the  notation  ru  -  dr/du,  rv  =  dr/dv  to  express  the  unit  tangent  vector  t  as 


where 


/  =  txru  +  t2rv , 


tx  _  du 
t2  dv 


Taking  the  dot  product  of  (22)  with  #*„  and  rv,  respectively,  we  obtain 


E  Fit, 


where 


F  G  t. 


E  =  ru-ru,  F  =  ru-rv,  G  =  rv  rv. 
From  (23)  and  (24),  we  then  obtain 

du  =  G(f  rv) 

dv  E(t  ■?')- Fit  -  F  ) 


If  this  yields  du/dv  <  1,  then  the  reciprocal  dv/du  should  be  used  to  initiate  the  iteration  process 
with  (17)  instead  of  (16).  Regardless,  knowledge  of  the  ratio  from  (26)  that  is  consistent  with 
movement  along  the  arc  in  the  ±  t0  direction  allows  setting  of  Auq  and  Avo  to  achieve  the  desired 
step  size,  \Aro\  =  Aso .  The  incremental  step  in  3-D  target  space  is  given  by 

A  r  =  ruAu  +  rvAv ,  (27) 

which  yields, 

(As)2  =  Ar  ■  Ar  =  E(Au)2  +  IFAuAv  +  G(Av)2 .  (28) 

Next,  setting  Au/Av  =  du/dv  or  Av/Au  —  dv/du  results  in 


Av  =  ±- 


+2F— +  G 
dv  J  dv 


(29a) 
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or 


A  u  —  ± 


As 


„  ^  „  dv  _ 

f^Vl 

,  E  +  2F —  +  G 

\(  du 

\^du  ) 

(29b) 


Finally,  select  the  sign  for  Av  (or  Au)  to  correspond  to  the  sign  for  t2  (or  ti)  obtained  from  (24). 
As  the  iteration  progresses,  arc  step  size  Asj  will  stray  from  the  original  and  desired  Aso-  This  is 
corrected  by  rescaling  Auj  and  Av,-  after  each  iteration  to  maintain  Aso-  Smaller  As  improves  the 
accuracy  of  the  numerical  solution  to  the  differential  equation. 


The  Christoffel  symbol  coefficients  [10,11]  in  (16)  and  (17)  follow: 


,  _GE,-2FFU  +  FE,  ,  _  2EF.  -  EE,  -  FE, 

11  2A  ’  11  2A 

(30a) 

i  _GEV-FGU  2_EGu-FEv 

12  2A  ’  12  2A 

(30b) 

,  2 GFV-GGU+FGV  2  _  EGV  - 2FFV  - FGU 

22  2A  ’  22  2A 

(30c) 

A  =  EG-F2, 

(30d) 

where  Eu  =  FE/du ,  etc. 


The  above  procedure  will  track  a  geodesic  across  a  single  patch  symbolically  and  parametrically 
specified  by  r  =  r(u,v).  When  the  geodesic  reaches  the  edge  of  a  patch,  it  is  necessary  to  continue 
the  geodesic  into  the  bordering  patch.  This  can  be  achieved  using  the  topology  map  in  a  manner 
similar  to  that  described  for  extending  the  SB  path  across  patch  borders. 

Figure  31  illustrates  the  implementation  of  the  creeping  wave  ray  tracer  on  a  cylinder  with 
hemisphere  end  caps.  The  creeping  rays  emanate  from  the  SB  in  the  direction  of  the  incident 
wave.  On  this  geometry,  one  can  observe  the  formation  of  both  caustic  points  and  the  more 
general  caustic  arcs  (lines).  This  issue  of  creeping  wave  caustics  is  discussed  below.  Note  also 
that  some  creeping  rays  emerge  back  into  the  lit  side  after  traveling  around  the  shadow  side.  For 
backscattering  applications,  these  rays  can  generate  significant  delayed  scattering  features  on 
smooth  bodies. 
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Incident 


Fig.  31.  Traced  creeping  wave  rays  on  a  cylinder  with  hemisphere  end  caps.  The  creeping  wave  ray  paths  are 
plotted  without  regard  to  energy  decay  along  the  geodesic  paths. 

The  next  section  will  discuss  the  means  of  assigning  proper  creeping  wave  surface  currents  to  the 
creeping  rays,  but  a  few  geometrical  points  remain  to  be  addressed.  The  creeping  rays  are  1-D 
arcs,  but  what  is  required  for  a  scattering  computation  are  creeping  ray  “tubes”  which  are  2-D 
ribbons  over  the  surface.  One  way  to  obtain  these  is  to  consider  two  creeping  ray  arcs  launched 
from  neighboring  points  on  the  SB  separated  by  a  sub-wavelength  distance  (e.g.,  A,/ 10).  The  two 
arcs  then  define  the  lateral  boundaries  of  a  creeping  wave  ray  tube.  Depending  on  local  surface 
curvature  conditions  encountered  as  these  geodesic  arcs  are  constructed,  the  ray  tube  walls  will 
either  converge  or  diverge.  When  traveling  over  flat  surfaces,  or  even  singly  curved  surfaces 
(e.g.,  cylindrical  surfaces),  the  divergence  rate  remains  fixed.  On  general,  doubly  curved  surfaces 
(e.g.,  sphere,  ellipsoid,  etc.),  divergence  rate  varies  as  the  ray  tube  propagates.  These  effects  can 
also  be  observed  in  Fig.  31. 

One  of  the  critical  geometrical  properties  required  in  evaluating  the  surface  currents  is  the 
divergence  factor  of  the  creeping  wave  ray  tube.  This  can  be  obtained  directly  by  measuring  the 
lateral  separation  between  the  creeping  ray  walls  at  any  point  along  their  respective  arcs.  This 
separation  must  be  determined  at  points  of  common  arc  length  (measured  from  the  SB)  between 
the  two  geodesic  arcs.  The  divergence  factor  is  given  by 


where  Wt  is  the  ray  tube  lateral  width  at  the  creeping  ray  sample  point  and  Wo  is  the  original 
width  at  the  SB. 
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Two  other  geometrical  quantities  are  required  by  the  creeping  wave  formulation:  arc  curvature 
and  arc  torsion.  The  arc  curvature  at  any  point  along  the  geodesic  is  given  by  [10] 


_  et\  +  “/V2  + 
arc  ~  Et]  +  2Ftf  +  Gt\  ’ 


(32) 


where  tl  and  t2  are  obtained  from  (24).  Since  tl  and  t2  are  also  proportional  to  Au  and  Av,  these 
can  be  substituted  directly  into  (32)  to  yield  Karc  as  the  geodesics  are  constructed.  The 
coefficients  E,  F  and  G  are  given  in  (25).  The  coefficients  e,  f,  and  g  are  given  by 


e  =  ruu-n,  f  =  rm-h,  g  =  rvv-h,  (33) 

where  h  is  the  2-D  surface  normal  and  ruu  =  3 r/du2,  etc. 

The  formulation  also  requires  the  arc  torsion  at  the  SB  launch  point  of  the  creeping  ray.  This  is 
given  by 


(34) 


where  K\  and  K2  are  the  principal  curvatures  of  the  surface.  They  are  the  two  roots  of  the 
quadratic  equation, 

(35) 


K2-2KmK  +  Kg=0, 


where 


_  eG-2/F  +  gE 
m  2(£G-F2)  ’ 


(36a) 


K  =  —  f  2. 
s  EG-F 2 


(36b) 


Km  and  Kg  are  termed  the  mean  surface  curvature  and  Gaussian  surface  curvature,  respectively. 


D.  Method  for  Determining  Creeping  Wave  Surface  Currents 

Direct  asymptotic  solutions  for  surface  currents  in  the  shadow  region  (and  the  transition  region 
around  the  SB)  are  not  available.  However,  researchers  [12,  13]  have  developed  workable 
formulations  for  arbitrary  convex  surfaces  by  heuristically  combining  asymptotic  solutions  to  the 
canonical  problems  of  metallic  sphere  and  cylinder  scattering  [2,14].  The  sphere  is  a  doubly 
curved  surface  with  creeping  rays  that  converge  toward  a  caustic  in  the  deep  shadow  region  and 
then  diverge  again.  The  cylinder  is  a  singly  curved  surface  with  non-diverging  creeping  waves 
that  also  introduces  the  effect  of  torsion  at  oblique  incidence  angles.  Hence,  these  two  canonical 
solutions,  together,  incorporate  the  key  physical  effects  for  a  heuristically  constructed  general 
solution. 
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1.  Fock  Creeping  Wave  Solution  for  a  Sphere 

To  anticipate  the  benefit  of  this  approach,  we  studied  scattering  based  on  the  creeping  wave 
formulation  for  the  PEC  sphere  (radius  a),  one  of  the  canonical  components  in  the  generalized 
heuristic  formulation.  For  this  case,  we  assume  a  plane  wave  incident  in  the  -z  direction  with  a 
unit  electric  field  polarized  in  the  +x  direction,  as  shown  in  Fig.  32. 


Fig.  32.  Incident  plane  wave  and  polarization  for  the  sphere  scatterer.  The  incident  direction  is  directly  toward  the 
sphere’s  north  pole,  and  the  shadow  boundary  will  be  the  equator  defined  along  the  xy -plane. 


The  total  magnetic  evaluated  at  the  sphere  surface  (r  =  a)  is  given  by  [2,13,14] 
H  =  He9  +  H/,  with 


H.= 


m 


and 


^ )]  -  yfe )] , 


where  0  <  9<  n  and 


D  = 


•s/sin  ~9  ’ 


/  1  ,n  1 2<  9  <n  (shadow  region) 
c0  =  1  , 

u / D  ,0  <9  < n/2 (lit region) 


m  = 


ka 


,1/3 


f  a{9  -n IT)  , shadow  region 

Tl  =  1  .  , 

-acos#-£,  13k  ,  lit  region 

|  mr,/fl  , shadow  region 
1  -mcos9  ,  lit  region 


(37a) 

(37b) 

(38a) 

(38b) 

(38c) 

(38d) 

(38e) 
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T2=a{37T/2-0),  (38f) 

%2=mr2/a.  (38g) 

In  these  expressions,  time-harmonic  dependence  exp(-jcot)  is  assumed  throughout.  Fock 
functions  f(£)  and  g(E)  are  given  by  [2,14]  as 


/fe)  = 


(39a) 

(39b) 


Where  wj(t)  and  its  derivative  are  Fock-type  Airy  functions,  r is  a  complex  contour  that  starts  at 
infinity  in  the  angular  sector  jt/3  <  arg  t  <  n  and  ends  at  infinity  in  the  angular  sector -jt/3  <  arg 
t  <  7tJ3.  Along  the  way,  /"passes  between  the  origin  and  the  pole  of  the  integrand  nearest  the 
origin.  The  properties  of f(E)  and  g(E)  and  methods  of  their  evaluation  are  available  from  [2,14], 

The  surface  current  can  be  obtained  from  (37)  using  J  =  n  x  H .  Alternatively,  using  the  PO 
approximation,  the  surface  currents  are  given  by  JP0  =2 nxH‘ .  The  resulting  surface  currents 
for  koa  -  3  are  shown  in  Fig.  33.  It  is  evident  that  the  creeping  wave  formulation  yields  more 
realistic  surface  currents  than  PO,  with  the  exception  of  the  singularity  at  0  =  n  (the  caustic), 
where  the  creeping  wave  formulation  is  known  to  break  down.  The  creeping  wave  formulation 
also  diverges  at  6=0  due  to  the  caustic  for  the  creeping  wave  traveling  around  the  back  side  and 
re-emerging  into  the  lit  region,  although  the  effect  is  less  pronounced  since  the  underlying 
creeping  wave  has  traveled  twice  as  far  as  the  distance  to  the  first  caustic. 
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Sphere  Surface  Currents:  kQa  -  3,  phi  -  0  deg 


Fig.  33.  Surface  currents  on  sphere  predicted  by  physical  optics  and  creeping  wave  formulation:  (a)  \Jg\  at  (/>=().  (b) 
\Jj  at  <j>=n/2.  Note  that  the  influence  of  the  caustic  points  at  the  sphere’s  poles  has  not  been  removed  in  these  plots. 


The  back-scattered  field  radiated  by  the  surface  currents  is  given  in  Fig.  34  as  a  frequency 
response.  This  is  the  same  plot  as  Fig.  27,  but  this  time  we  include  the  result  for  the  creeping 
wave  formulation.  Clearly,  the  creeping  wave  formulation  is  a  significant  improvement  over  PO, 
using  the  exact  solution  Mie  series  as  the  reference. 

The  time  domain  (range  profile)  response  of  the  sphere  (< a  =  0.0477  m )  in  Fig.  35  also  shows 
significant  improvement  using  the  creeping  wave  formulation.  The  false  scattering  feature 
predicted  by  PO  at  down  range  position  z  =  0  has  been  completely  eliminated  owing  to  the 
uniform  expressions  for  the  surface  currents  current  in  the  vicinity  of  the  SB  in  (37).  Also,  the 
delayed  response  at  z  =  m/2  =  0.075  m,  due  to  energy  traveling  around  the  back  side  of  the 
sphere  and  re-emerging  on  the  other  side,  is  also  correctly  predicted.  However,  the  creeping 
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wave  formulation  also  introduces  a  false  scattering  feature  at  z  =  a(n/2  +  l)/2  =  0.061  m.  This  is 
caused  by  the  breakdown  of  the  formulation  in  the  vicinity  of  the  caustic  at  0  =  n. 


Frequency  Response  of  PEC  Sphere 


Fig.  34.  Frequency  back-scattering  response  of  the  PEC  sphere  computed  by  exact  solution  (Mie  series),  asymptotic 
creeping  wave  formulation,  and  physical  optics. 


Range  Profile  of  PEC  Sphere  (a  -  0.0477  m) 


Fig.  35.  Time-domain  back-scattering  response  of  PEC  sphere  ( a  =  0.0477  m)  computed  by  exact  solution  (Mie 
series),  asymptotic  creeping  wave  formulation,  and  physical  optics. 


To  keep  the  surface  current  radiation  integral  from  diverging,  the  integration  is  abruptly 
terminated  at  0=1 79°.  This  gives  a  better  result  than  integrating  all  the  way  through  caustic  at  0 
=  180°,  but  a  false  scattering  feature  is  still  generated  by  the  abrupt  truncation.  There  are 
convergent  asymptotic  expressions  for  the  surface  currents  in  the  vicinity  of  the  caustic  on  the 
sphere  [2,14],  but  these  do  not  so  readily  lend  themselves  to  a  generalized  solution  for  arbitrary 
curved  surfaces.  In  practical  problems  involving  targets  with  doubly  curved  surfaces,  caustics 
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will  occur,  so  this  problem  must  be  confronted.  A  practical  strategy  may  be  to  smooth  the 
currents  around  the  caustic  so  as  to  avoid  both  the  caustic  singularity  and  the  false  scattering 
features  created  by  abrupt  truncation  of  the  radiation  integral.  At  the  time  of  this  reporting, 
however,  this  approach  has  not  been  implemented  or  tested  in  our  study. 

2.  General  Solution  for  Creeping  Waves  in  the  Shadow  Region 

For  a  general  creeping  wave  expressions  that  can  be  used  with  the  arbitrary  geodesic  paths  over 
convex  surfaces,  it  is  useful  to  define  a  local  right-handed  coordinate  system,  n,b,t ,  for  the 
creeping  rays,  as  shown  in  Fig.  36.  At  all  points  along  the  geodesic  arc,  h  is  the  outward  surface 
normal,  and  unit  tangent  vector  t  is  in  the  direction  of  creeping  ray  travel.  The  bi-normal  vector 
is  given  by  b  =  t  x  h ;  it  tracks  the  torsion  of  the  geodesic.  The  expressions  for  the  surface  fields 

/s  A 

and  surface  currents  created  by  the  creeping  waves  will  be  given  in  terms  of  the  local  b  and  t 
unit  vectors. 


Fig.  36.  Local  coordinate  system  and  notation  for  geodesic  creeping  ray  arcs. 


The  creeping  ray  launches  from  the  SB  at  point  Q\  also  called  the  source  point.  Arbitrary  points 
along  the  geodesic  arc  are  labeled  Q.  In  expressions  given  below,  quantities  like  K(Q)  mean  that 
function  evaluated  at  point  Q  on  the  geodesic  creeping  ray  arc. 

Without  repeating  the  derivation  in  Mitschang  and  Wang  [13]1,  the  formulas  for  the  creeping 
wave  surface  magnetic  fields  in  the  shadow  region  are  given  by 


Hsh=Hbb  +  H,t 


(40a) 


Hb=(n-ei 


1  ¥,(Q) 


(40b) 


1  Mitchang  and  Wang  [13]  assume  exp(jcot)  time-harmonic  dependence  in  their  expressions,  although  their  paper 
mistakenly  claims  exp  (-jot)  dependence  at  one  point.  This  report  assumes  exp(-jat)  dependence,  so  all  phasor 
expressions  appear  as  the  complex  conjugate  of  their  correspondents  in  [13]. 
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f2 +j 


A(g)a;(e) 

kps(Q) 


(40c) 


where  e‘  and  h‘  are  the  incident  plane  wave  polarization  unit  vectors  for  the  electric  and  magnetic 
fields,  respectively.  Parameter  A  is  the  transition  parameter  that  heuristically  combines  the 
canonical  sphere  and  cylinder  asymptotic  creeping  wave  solutions  to  yield  a  more  general 
solution  for  arbitrary  convex  surfaces,  A  is  given  by 


(  0  ,  cylinder 

=  1  ,  sphere 

R2(Q')/ (Q'),  general  convex  surfaces 


(41) 


where  Ri(Q’)  and  Ri(Q’)  are  the  principle  radii  of  curvature  at  the  SB  launch  point  of  the 
creeping  ray;  Ri(Q’)  >=  RjiQ’).  Note  that  A(Q’)  is  a  fixed  quantity  based  on  surface  conditions 
at  the  SB,  not  a  varying  quantity  along  the  creeping  ray. 

Factor  D(Q)  in  (40)  is  the  divergence  factor  of  the  creeping  ray  tube.  A  method  for  obtaining 
D(Q)  at  any  point  along  the  creeping  ray  was  discussed  earlier.  In  general,  D(Q)  is  a  complex 
number,  as  discussed  further  in  the  section  below  on  caustic  crossings.  pg(Q)  is  the  geodesic 
radius  of  curvature  at  any  point  along  the  geodesic  arc.  It  is  obtained  directly  from  (32)  since 
Pg(Q)  =  l/Karc(Q).  Factor  T0(Q’)  in  (40b)  is  directly  related  to  the  creeping  ray  torsion  T(Q’)  at 
the  SB:  Tq(Q’)  =  T(Q’)  pg(Q’).  T(Q’)  is  obtained  from  (34). 


Factors  Fj  -  Fj(Q,Q’)  and  F2  =  F2(Q,Q’)  incorporate  the  Fock  functions  f(S)  and  g(£)  governing 
the  decay  of  creeping  waves  as  they  extend  from  the  shadow  boundary.  Fj  and  F2  are  given  by 


F,(Q,Q')=e,k’ 


M 


nl/6 


pM 


d(q fete). 


AQ.Q!)=^-e 


-  I  Jkt 

m(Q) 


pM 

p,(Q ) 


— il  /  6 


(42a) 


(42b) 


t-  v(Q,Q’)  is  the  arc  length  of  the  creeping  ray  measured  from  the  SB.  Expressions  for  the  Fock 
functions  f(S)  and  g( E)  are  given  in  (39).  Fock  parameter  %=  %(Q,Q’)  is  given  by 


(43) 
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where 


m(s)  =  m{Q)  = 


kpMT 


(44) 


3.  Caustic  Crossing 

For  the  sphere  geometry,  the  creeping  waves  reach  a  caustic  at  0  =  n.  This  is  the  point  where  the 
creeping  waves  launched  at  the  SB  (at  6  =  n/2)  converge  at  the  pole  of  the  sphere.  Beyond  the 
caustic,  the  creeping  waves  will  again  diverge  until  they  reach  the  SB  on  the  other  side  of  the 
sphere  from  which  they  were  launched.  Beyond  that  point,  they  will  again  converge  to  a  caustic, 
this  time  at  6=0.  On  the  sphere,  this  process  continues  indefinitely  as  the  creeping  waves 
continue,  forming  successive  caustics  at  6  =  ;rand  6=0. 


In  the  general  implementation  of  creeping  waves  on  doubly  curved  surfaces,  caustics  will  arise. 
In  general,  they  do  not  form  points.  Rather,  neighboring  creeping  waves  rays  will  converge  over 
a  locus  of  points  forming  a  curved  line  on  the  surface.  This  caustic  line  is  the  1-D  correspondent 
to  the  2-D  caustic  surfaces  that  arise  in  geometrical  optics  (GO)  with  the  reflection  of  rays  off  of 
curved  surfaces  [15].  As  in  the  GO  case,  the  fields  predicted  by  the  asymptotic  creeping  wave 
solution  are  invalid  in  the  immediate  vicinity  of  the  caustic,  but  fields  both  before  and  beyond  the 
caustic  are  valid.  Also,  in  crossing  the  caustic,  a  phase  shift  occurs  that  needs  to  be  accounted  for 
in  the  numeric  implementation. 


If  one  inspects  the  canonical  solution  of  the  sphere  in  (37),  which  is  valid  for  0  <  6  <  n,  the 
expression  includes  extra  terms  for  creeping  waves  that  start  at  the  opposite  SB  and  travel  around 
the  back;  these  are  the  terms  with  subscript,  2.  This  is  in  contrast  to  (40),  which  does  not 
implicitly  account  for  creeping  waves  arriving  from  a  different  SB  launch  point.  Doing  so  is  not 
possible  on  general  convex  surfaces,  and  the  intent  with  the  general  implementation  is  for  this 
effect  to  be  created  explicitly  by  launching  creeping  wave  rays  from  all  points  along  the  SB  and 
allowing  them  to  travel  freely. 

Returning  to  the  sphere  expressions  in  (37),  and  casting  these  in  terms  of  the  general  solution  of 

/v  A 

(40)  and  its  n,b,t  framework  for  specifying  field  and  surface  current  directions,  the  phase  shift 
for  a  creeping  wave  crossing  a  caustic  can  be  deduced.  Accounting  for  the  flip  in  the  polarization 
definitions  of  6  and  <f)  at  the  poles,  and  also  the  relative  change  in  orientation  of  the  incident 
fields  at  the  opposite  SB  point,  it  can  be  shown  from  (37)  that  the  terms  proportional  to 
divergence  factor  D  experience  a  -9(f  (i.e.,  - j )  phase  shift  when  crossing  the  caustic,  while  terms 
proportional  to  D3  experience  a  +90°  (i.e.,  +j)  phase  shift. 


As  with  the  GO  caustics,  the  creeping  wave  caustic  phase  shifts  can  be  interpreted  as  arising 
from  a  generalized  divergence  factor  D  that  takes  on  complex  values.  In  (38a),  D  is  real  and 
positive  because  0  <  6  <  n.  If  we  recast  the  backside  creeping  wave  terms  in  (37)  that  travel 
through  the  caustic  to  use  the  range,  n  <  6  <  2n,  in  their  various  arguments,  then  D  becomes 
imaginary  and  (38a)  can  be  re-written  as, 
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(45) 


Vsin#  has  branch  points  at  6  =  nn,  n  =  {...,-2, -1,0,1, ...}. 

If  we  choose  the  branch  where 
V sin 9  =  /'Vsin#' ,  0=  0’  +  nn,  0  <  6’  <  n, 

then  both  D  and  D3  take  on  the  correct  real  or  imaginary  sign,  {+l,+j,-l,-j},  in  the  generalized 
version  of  (37).  For  the  generalized  implementation  of  Eqs.  (40)  -  (44),  whenever  a  caustic  is 
encountered  in  evaluating  D  in  (31),  a  -j  factor  should  be  applied  to  D  while  traveling  forward 
through  the  caustic.  For  the  D3  terms,  this  factor  should  be  included  before  cubing  D. 

4.  Uniform  Solution  For  Lit  Region 

Equation  (40)  describes  the  creeping  wave  surface  currents  in  the  shadow  region.  If  these  are 
combined  with  the  PO  surface  currents  in  the  lit  region,  there  will  be  an  abrupt  discontinuity  at 
the  SB.  This  is  because  PO  breaks  down  near  the  SB.  Since  one  of  the  important  goals  of  the 
creeping  wave  formulation  is  to  improve  over  PO  near  the  SB,  not  just  in  the  deep  shadow 
region,  researchers  have  developed  generalized  expressions  based  on  the  Fock’s  asymptotic 
solution  that  can  be  used  in  the  lit  region  and  create  a  smooth  transition  over  the  SB.  For  the  lit 
region  coordinate  system  shown  in  Fig.  37,  Mitschang  and  Wang  [13]  give  the  following 
expression  for  fields  in  the  lit  region: 

Hlit  =  {e‘  ■  n]F;“b  -  (h‘  ■  h^t ,  (46) 

where 

if=e'*"gte),  (47a) 

(47b) 

m,\Q) 

and 


(48a) 

(48b) 

(48c) 
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Note  that  in  Equation  (48b),  the  torsion-related  parameter  is  evaluated  at  the  point  of  incidence  Q 
rather  at  the  shadow  boundary  Q  This  is  because  we  are  in  the  lit  region  and  all  points  are 
“source”  points.  Also,  the  quantity  m(Q)  in  (48b)  is  evaluated  using  (44)  based  on  the  curvature 
along  the  tangent  direction  defined  in  Fig.  37. 

It  should  be  pointed  out  that  the  expression  in  (46)  lacks  the  second-order  (1/k)  terms  present  in 
(40)  arising  from  the  sphere  solution  in  (37).  Hence,  there  will  still  be  a  small  discontinuity  in 
the  surface  currents.  This  effect  is  evident  in  the  numerical  results  we  present  later  for  the  sphere 
using  the  generalized  formulation.  Also,  Pathak  [12]  gives  an  expression  that  they  claim 
maintains  continuity  at  the  SB,  and  which  includes  additional  terms  relating  to  the  surface  torsion 
beyond  those  given  in  (46).  At  the  time  of  this  reporting,  we  are  conducting  further  investigation 
into  this  issue  to  resolve  these  discrepancies  and  obtain  a  uniform  solution  across  the  SB. 

E.  Results  and  Discussion 

To  test  the  reasonableness  of  the  formulation  for  the  shadow  region  creeping  wave  surface 
currents,  the  study  code  was  employed  to  predict  results  on  a  sphere  and  cylinder  with 
hemisphere  end  caps.  The  surface  currents  are  shown  on  the  target  bodies  in  Figs.  38  and  39. 
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Fig.  38.  Predicted  shadow  region  creeping  wave  currents  on  a  sphere.  Only  a  subset  of  the  creeping  waves  are 
shown  launching  from  a  roughly  170°  portion  of  the  incident  shadow  boundary  in  order  to  clearly  illustrate  currents 
in  both  the  near  and  far  shadow  regions. 


Fig.  39.  Predicted  shadow  region  creeping  wave  currents  on  a  cylinder  with  hemisphere  caps.  For  clarity,  only  a 
subset  section  of  the  creeping  wave  are  shown. 


Inspecting  the  surface  currents,  several  things  are  evident.  First,  the  magnitudes  of  the  currents 
generally  decrease  as  they  travel,  corresponding  to  their  progressive  radiation  to  space. 
Mathematically,  this  arises  directly  from  the  large  positive  argument  behavior  of  the  Fock 
functions.  Also,  as  the  surface  currents  approach  a  caustic,  their  intensity  increases  toward  a 
singularity;  this  is  most  evident  for  the  sphere  geometry  in  Fig.  38.  The  singularity  arises  directly 
from  the  divergence  factor.  Surface  currents  immediately  near  the  caustics  are  not  shown  (i.e., 
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painted  black)  in  Figs.  38  and  39.  This  corresponds  directly  to  the  logic  in  the  study  code  for 
skipping  the  surface  integration  near  the  caustic  to  avoid  a  divergent  result.  In  the  case  of  the 
sphere  in  Fig.  38,  the  geodesics  converge  on  a  point  caustic.  That  is  a  special  case.  For  the  more 
general  shape  in  Fig.  39,  the  loci  of  the  caustics  are  curved  lines.  Finally,  the  representation  of 
surface  current  phase  reveals  the  expected  steady  progression  arising  directly  from  the 
exponential  factors  in  (42). 

To  further  verify  the  numeric  implementation  of  the  general  formulation,  scattering  computations 
were  performed  on  the  sphere  geometry  (radius  a  =  2  m)  and  compared  to  both  the  exact  Mie 
series  solution  and  physical  optics.  The  time  response  (i.e.,  range  profile)  is  given  in  Fig.  40  and 
the  frequency  response  in  Fig.  41 


Range  Profile:  PEC  Sphere  (a  -  2  m) 


Fig.  40.  Back-scattering  time  response  (range  profile)  for  a  PEC  sphere  (radius  a  =  2  m)  computed  from  exact 
solution  (Mie)  series,  physical  optics,  and  numeric  implementation  of  the  generalized  creeping  wave  ray  tracer. 
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Fig.  41.  Back-scattering  frequency  response  for  a  PEC  sphere  (radius  a  =  2  m)  by  exact  solution,  physical  optics,  and 
generalized  creeping  wave  ray  tracer. 

For  the  range  profile  in  Fig.  40,  we  observe  that  the  generalized  implementation  reintroduces  the 
false  scattering  feature  at  z  =  0,  as  compared  with  the  more  accurate  results  based  on  the 
specialized  formulation  for  the  sphere  (37)  in  Fig.  35.  This  is  caused  by  the  remaining 
discontinuity  in  surface  currents  around  the  SB  arising  from  inconsistencies  between  the  shadow 
region  formation  of  (40)  and  the  lit  region  formulation  of  (46).  This  was  discussed  earlier,  and 
our  current  work  centers  on  reconciling  these  two  representations.  Also,  the  numeric 
implementation  of  (46)  in  the  study  code  has  not  been  firmly  established  by  our  own  inspection 
of  the  source  code  at  this  reporting.  Nonetheless,  even  in  the  current  implementation,  this 
spurious  feature  is  reduced  relative  to  PO. 

Because  of  the  limited  bandwidth  in  the  computed  frequency  response,  in  the  range  profile  from 
the  general  creeping  wave  implementation,  the  false  scattering  feature  at  z  =  a(n/2  +  l)/2  =  2.57 
m  is  not  completely  resolved  from  the  trailing  response  at  z  =  an/2  =  3.14  m  due  to  emergence 
creeping  wave  back  into  the  lit  region.  However,  its  effect  is  observed  in  the  widened  and 
distorted  response  feature  in  the  vicinity  of  these  two  locations.  This  spurious  feature  was  also 
evident  in  the  higher-resolution  range  profile  based  on  the  specialized  formulation  of  the  sphere 
in  Fig.  35.  As  discussed  earlier,  this  is  caused  by  abruptly  halting  the  surface  integration  in  the 
vicinity  of  the  caustic  at  the  pole.  This  situation  should  be  improved  by  instead  smoothing  the 
surface  currents  around  the  caustic  to  avoid  a  non-integrable  singularity  (caused  by  breakdown  of 
the  creeping  wave  formulation),  but  this  remains  to  be  tested  numerically.  Finally,  the  trailing 
response  due  to  the  emergence  of  the  creeping  waves  back  into  the  lit  side  after  traveling  around 
the  back  is  properly  produced  by  the  generalized  implementation  of  the  asymptotic  creeping 
wave  formulation. 

Correspondingly,  the  frequency  response  in  Fig.  41  also  suffers  in  the  generalized 
implementation.  This  is  caused  by  the  same  effects  just  discussed.  In  spite  of  these  remaining 
problems,  it  still  represents  an  improvement  over  the  PO  solution  in  terms  of  picking  up  the 
higher  rate  of  oscillation  about  the  final  RCS  level  converged  upon  for  k0a  »  1.  If  the 
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remaining  issues  of  the  SB  discontinuity  can  be  resolved,  the  generalized  creeping  wave 
implementation  holds  strong  promise  for  greatly  improving  upon  the  base  PO  solution  of 
Xpatch®. 

F.  Summary 

A  complete  method  was  presented  for  improving  the  performance  of  Xpatch®  in  the  deep  shadow 
region  on  smoothly  curved  bodies  using  a  creeping  wave  approach  that  is  consistent  with  the 
Xpatch  framework  of  approximating  surface  currents  on  arbitrary  target  geometries  through  ray 
tracing.  This  included  an  alternate  technique  for  identifying  the  shadow  boundary  on  the  target, 
which  was  qualitatively  compared  with  the  original  SB-finding  technique  developed  in  this 
project  for  implementation  of  the  SB-ELDCs.  The  overall  method  also  included  numeric 
approaches  for  tracing  the  geodesic  creeping  wave  rays  from  the  shadow  boundary  and 
developing  appropriate  surface  currents  along  them  from  asymptotic  formulations  based  on  the 
Fock  theory. 

These  methods  were  implemented  in  a  research  code  employed  by  this  study  and  tested  on  simple 
cases.  Based  on  this  evaluation,  we  reach  the  following  conclusions.  First,  the  alternate  method 
for  finding  the  SB  appears  workable,  and  is  well-suited  toward  the  application  of  launching 
creeping  waves  in  the  shadow  region  since  it  constructs  the  SB  as  a  contiguous  curve  that  can  be 
re-sampled.  However,  further  work  is  required  in  adapting  it  to  more  complex  geometries  that 
contain  multiple  SBs.  Second,  the  method  for  tracing  the  geodesic  creeping  wave  rays  behaved 
both  reasonably  and  efficiently  on  the  curved  surface  geometries  considered.  This  included  the 
ability  to  carry  the  geodesics  across  adjoining  patches  whose  topological  links  were  not 
established  a  priori  as  input  to  the  code,  an  important  practical  benefit  to  the  method  described. 
Third,  and  finally,  the  method  for  applying  shadow  region  surface  currents  to  the  creeping  wave 
rays  improved  the  scattering  solution  on  the  sphere  test  case  considered.  However, 
discontinuities  in  the  lit  side  and  shadow  side  formulations  remain  that  still  cause  undesirable 
scattering  features  in  the  generalized  implementation.  Our  current  work  focuses  on  examining 
this  issue  further  so  that  these  non-physical  SB  discontinuities  in  the  surface  currents  can  be 
eliminated. 
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VI.  CONCLUSIONS 

The  Xpatch®  code  with  the  SB-ILDCs  shows  a  large  improvement  in  calculating  the  scattering 
from  a  number  of  objects.  Large  cylinders  and  spheres  were  modeled.  The  effectiveness  of  the 
SB-ILDCs  was  illustrated  for  objects  with  different  sized  radii  of  curvature  and  for  objects  that 
possess  a  varying  radii  of  curvature. 

The  Xpatch®  code  with  the  TW-ILDCs  calculates  with  greater  accuracy  the  scattering  from 
objects  with  edges.  The  new  TW-ILDCs  eliminate  the  integration  singularities  present  in  the  old 
formulation  and  produce  scattering  results  that  are  free  of  the  previous  non-physical  field 
discontinuities.  The  original  Xpatch®  formulation  lacked  cross-polarization  terms  (infinity 
avoidance).  The  TW-ILDCs  are  based  on  a  rigorous  derivation  that  includes  cross  terms.  The 
TW-ILDCs  also  achieve  a  partial  accounting  of  2nd  order  diffraction  effects. 

A  detailed  study  of  deep  shadow  creeping  wave  scattering  as  a  means  of  supplementing  the  SB- 
ILDC  work  was  performed.  This  study  employed  a  prototype  algorithm  for  tracing  creeping 
waves  along  smoothly  curved  bodies.  This  algorithm  consists  of  a  geometrical  component  in 
which  contiguous  shadow  boundary  contours  are  extracted  followed  by  geodesic  ray  tracing 
component  to  determine  the  paths  of  creeping  waves  launched  at  the  shadow  boundary.  Details 
of  the  underlying  creeping  wave  EM  radiation  formulation  were  presented  along  with  numerical 
results  for  the  sphere  and  hemisphere-capped  circular  cylinder.  Comparisons  with  closed  form 
and  full  wave  results  are  good  and  demonstrate  the  added  physics  due  to  the  late  time  scattering 
of  creeping  waves.  However,  the  underlying  EM  approach  possesses  significant  shortcomings 
with  respect  to  caustic  points  along  the  geodesic  path.  Future  research  will  explore  development 
within  an  ILDC  framework  in  which  direct  calculation  of  the  surface  currents  can  be  avoided  in 
favor  of  more  stable  line  integral  representations  of  the  scattering. 
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